! This is a 1D advection example using square initial condition and periodic
! boundary condition for FFSL finite volume scheme.
!
! Li Dong <dongli@lasg.iap.ac.cn>
!
! - 2018-03-22: Initial creation.
! - 2024-03-22: Use adv_1d_square_case_mod.
! - 2024-03-24: Put ppm and slope into ffsl_common_mod.

program ffsl_adv_1d_case

  use adv_1d_test_case_mod
  use ffsl_common_mod

  implicit none

  real, allocatable :: rho(:,:)           ! Tracer density being advected at cell centers
  real, allocatable :: flx(:)             ! Flux at cell interfaces
  integer, parameter :: ns = 3            ! Stencil width
  integer i
  character(30), parameter :: scheme = 'ffsl_1d'

  namelist /params/ nx, nt, dt, flux_type, limiter_type, u

  call get_command_argument(1, namelist_path)
  inquire(file=namelist_path, exist=is_exist)
  if (is_exist) then
    open(10, file=namelist_path)
    read(10, nml=params)
    close(10)
  end if

  allocate(rho(1-ns:nx+ns,2))
  allocate(flx(1-ns:nx+ns))

  call adv_1d_test_case_init('square', ns, rho(:,old))
  call output(scheme, time_step, ns, nx, x, rho(:,old))

  ! Run integration.
  write(*, *) time_step, sum(rho(1:nx,old))
  do while (time_step < nt)
    call ffsl(rho(:,old), flx)
    do i = 1, nx
      rho(i,new) = rho(i,old) - (flx(i) - flx(i-1))
    end do
    call apply_bc(ns, nx, rho(:,new))
    call advance_time()
    call output(scheme, time_step, ns, nx, x, rho(:,old))
    write(*, *) time_step, sum(rho(1:nx,old))
  end do

  call adv_1d_test_case_final()

  deallocate(rho)
  deallocate(flx)

contains

  subroutine ffsl(q, f)

    real, intent(in ) :: q(1-ns:nx+ns)
    real, intent(out) :: f(1-ns:nx+ns)

    real cfl, cf, ql, dq, q6
    real s1, s2, ds1, ds2, ds3
    integer i, j, ci, iu

    !  Cell:  |--0--|--1--|--2--|--3--|
    !  Face: -1     0     1     2     3
    face_loop: do i = 0, nx
      f(i) = 0
      cfl = u * dt / dx
      ci = int(cfl)
      cf = cfl - ci
      ! Calculate integer flux.
      if (cfl > 0) then
        do j = 1, ci
          f(i) = f(i) + q(i - j + 1)
        end do
      else if (cfl < 0) then
        do j = 1, -ci
          f(i) = f(i) - q(i + j)
        end do
      end if
      iu = merge(i - ci, i - ci + 1, cfl > 0)
      ! Calculate fractional flux.
      select case (flux_type)
      case ('upwind')
        f(i) = f(i) + cf * q(iu)
      case ('van_leer')
        dq = slope(q(iu-1), q(iu), q(iu+1))
        f(i) = f(i) + cf * (q(iu) + (sign(1.0, cfl) - cf) * dq * 0.5)
      case ('ppm')
        call ppm(q(iu-2), q(iu-1), q(iu), q(iu+1), q(iu+2), ql, dq, q6)
        if (cfl >= 0) then
          s1 = 1 - abs(cf)
          s2 = 1
        else if (cfl < 0) then
          s1 = 0
          s2 = abs(cf)
        end if
        ds1 = s2    - s1
        ds2 = s2**2 - s1**2
        ds3 = s2**3 - s1**3
        f(i) = f(i) + sign(1.0, cfl) * (ql * ds1 + 0.5 * dq * ds2 + q6 * (0.5 * ds2 - ds3 / 3.0))
      end select
    end do face_loop

  end subroutine ffsl

end program ffsl_adv_1d_case
